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ABSTRACT 

Using high resolution hydrodynamical simulations, we explore the spin evolution of 
massive dual black holes orbiting inside a circumnuclear disc, relic of a gas-rich galaxy 
merger. The black holes spiral inwards from initially eccentric co- or counter-rotating 
coplanar orbits relative to the disc's rotation, and accrete gas that is carrying a net 
angular momentum. As the black hole mass grows, its spin changes in strength and 
direction due to its gravito-magnetic coupling with the small-scale accretion disc. We 
find that the black hole spins loose memory of their initial orientation, as accretion 
torques suffice to align the spins with the angular momentum of their orbit on a 
short timescale (< 1 — 2 Myr). A residual off-set in the spin direction relative to the 
orbital angular momentum remains, at the level of < f 0° for the case of a cold disc, 
and < 30° for a warmer disc. Alignment in a cooler disc is more effective due to the 
higher coherence of the accretion flow near each black hole that reflects the large-scale 
coherence of the disc's rotation. If the massive black holes coalesce preserving the 
spin directions set after formation of a Keplerian binary, the relic black hole resulting 
from their coalescence receives a relatively small gravitational recoil. The distribution 
of recoil velocities inferred from a simulated sample of massive black hole binaries 
has median < 70kms _1 , much smaller than the median resulting from an isotropic 
distribution of spins. 

Key words: black hole physics - hydrodynamics - galaxies: starburst - galaxies: 
evolution - galaxies: nuclei 



1 INTRODUCTION 

The massive black holes (MBHs) that we observe today in lo- 
cal spheroids (Ferrarese & Ford 2005, and references therein) 
are expected to have grown through a series of major ac- 
cretion episodes in symbiosis with the growth of their host 
galaxies. Gas-rich major mergers may be at the heart of this 
joint evolution as they may explain the morphology of the 
hosts and at the same time account for the fueling of the un- 
derlying MBHs (e.g. Di Matteo, Springel & Hernquist 2005, 
and references therein). In the currently favored cold dark 
matter hierarchical cosmologies, galaxy mergers play indeed 
a key role in growing galaxies to their present sizes and the 
coalescence of MBHs in binaries is therefore expected to 
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be relatively common (Menou, Haiman & Narayanan 2001; 
Volonteri, Haardt & Madau 2003). 

Following a galaxy major merger, the pair of MBHs 
first interacts with stars (Begelman, Blandford & Rees 1980; 
Milosavljevic & Merritt 2001) and gas (Mayer et al. 2007). 
The pair loses orbital angular momentum, and it is ex- 
pected to harden progressively down to subparsec scales 
where emission of gravitational radiation drives the MBH 
inspiral down to coalescence. Depending on the properties of 
the coalescing binary, the pattern of the gravitational wave 
emission can be anisotropic, resulting in a non zero recoil 
velocity (a "kick" ) of the MBH remnant (Redmount & Rees 
1989). Several attempts to compute analytically the strength 
of the kick have been undertaken (Peres 1962; Bekenstein 
1973; Fitchett 1983; Fitchett & Detweiler 1984; Redmount 
& Rees 1989; Wiseman 1992; Favata, Hughes & Holz 2004; 



2 Dotti et al. 



Blanchet, Qusailah & Will 2005; Damour & Gopakumar 
2006; Schnittman & Buonanno 2007. 

Recent numerical simulations of the coalescence of spin- 
ning MBHs in full general relativity have been able to calcu- 
late explicitly kick velocities for a series of different binary 
configurations. It is found that three parameters influence 
the magnitude of the gravitational recoil of the relic MBH: 
the binary mass ratio, the spins, and the mutual orientation 
of the spins with respect to the orbital angular momentum. 
The recoil is largest, up to 4000 km s~ 1 , for nearly equal mass 
MBHs with large spins, when the spin vectors have oppo- 
site directions and are in the orbital plane (Campanelli et 
al. 2007). By contrast, recoils of < 200 km s -1 are imparted 
to the MBH remnant if the spins of the progenitors prior 
coalescence are orthogonal to their orbital plane. 

Purely general relativistic (GR) effects (i.e., spin-orbit 
and spin-spin interactions) may produce low recoil config- 
urations, when the MBH pairing is driven by gravitational 
wave emission. However, those GR effects depend strongly 
on the initial relative orientation of the MBH spins, and can 
result in low recoil configurations only for a small region of 
parameter space (Schnittman 2004; Herrmann et al. 2009; 
Lousto et al. 2009). 

In gas-rich mergers between galaxies of comparable 
mass (i.e. major mergers), close binary MBHs form under 
the action of dynamical friction against the gaseous and stel- 
lar background (Mayer et al. 2007; Callegari et al. 2009; see 
Colpi & Dotti 2009 for a review). During their inspiral MBHs 
are surrounded by a dense cocoon of gas that drives their dy- 
namical decay and provides fuel for the feeding of the holes 
(Dotti et al. 2007, 2009). Since matter carries angular mo- 
mentum also the spin vector can change during the accretion 
process. The details of the dynamics may have a profound 
influence on the mass and spin evolution of the two MBHs, 
and thus on the recoil velocity of the MBH resulting from 
their coalescence, and this is matter of our concern in this 
paper. 

The spin evolution during the MBH inspiral in a gas 
rich merger remnant has many implications . Spins, prior to 
coalescence, influence the extent of the gravitational recoil, 
and so the retention of the relic MBH inside its host galaxy. 
Accordingly, the spin distribution of the coalescing binaries, 
is critical as it determines the frequency of MBH retention in 
the host halo (Volonteri & Rees 2006; Volontcri 2007; Volon- 
teri, Haardt & Gultekin 2008; Gultekin et al. in prepara- 
tion). The magnitude of the MBH spins and their orientation 
relative to the orbit during mergers is also critical in shaping 
the stellar density profiles in ellipticals, as the kicked MBH 
moving on a return orbit can deposit its excess kinetic energy 
into the stellar background, causing the formation of stellar 
core (Boylan-Kolchin, Ma & Quataert 2004; Gualandris & 
Merritt 2008). A recoiling MBH can have an observational 
signature when moving across the host galaxy, creating an 
X-ray tail in the perturbed hot gas (Devecchi et al. 2009), 
an off-set active nucleus (Loeb 2007; Volonteri & Madau 
2008), shocking the inner rim of the accretion disc (Lippai, 
Frei & Haiman 2008; Schnittman & Krolik 2008), or drag- 
ging a stellar cusp with peculiarly high velocity dispersion 
(Merritt, Schnittman & Komossa 2009). Furthermore, spin 
orientations have important implications for gravitational 
wave astronomy, and for using gravitational wave measure- 
ments to constrain the formation history of MBHs (Vecchio 



2004; Lang & Hughes 2006; Berti & Volonteri 2008; Arun et 
al. 2009a, 2009b). 

Bogdanovic, Reynolds & Miller (2007) proposed a phys- 
ical process that could align the MBH spins with the orbital 
angular momentum of the binary, thus leading to slow re- 
coils for the MBH remnant. The key process for alignment 
is the presence of a coherent gas inflow. They speculate that 
accreting gas exerts gravito-magnetic torques that suffice to 
align the spins of both the MBHs with the angular momen- 
tum of the large-scale gas flow in which the orbit is embed- 
ded. 

Spin-disc alignment due to gravito-magnetic coupling 
has been studied by a number of authors in the case of iso- 
lated MBHs surrounded by their own discs (Bardeen & Pet- 
terson 1975; Natarajan & Pringle 1998; Scheuer & Feiler 
1996, Martin, Pringle & Tout 2007; Perego et al. 2009). Here 
we attempt to explore for the first time spin-disc alignment 
around MBH binaries. We expect low recoils when the spin- 
disc coupling is strong, i.e. when: 

• The two MBHs accrete > 1% of their initial mass before 
the coalescence (Natarajan & Pringle 1998; Natarajan & Ar- 
mitage 1999; Volonteri, Sikora & Lasota 2007; Perego et al. 
2009); 

• The accreted gas carries angular momentum in a preferred 
direction, flowing onto each MBH along a preferential plane 
determined by the distribution of angular momentum of the 
gas in the environment of the MBH. 

The first requirement can be fulfilled if the MBHs pair 
inside a dense, massive gaseous nuclear disc (Dotti et al. 
2009), such as that predicted to form in remnants of gas- 
rich major mergers by Mayer et al. (2007). To constrain the 
second requirement, we analyse a set of 3D Smooth Parti- 
cle Hydrodynamics (SPH) simulations already discussed in 
Dotti et al. (2009). The high resolution of these simulations 
enables us to resolve the gravitational sphere of influence 
of each MBH during their inspiral inside the circumnuclear 
disc, and to map the distribution of angular momentum of 
the SPH particles in the MBH vicinity. MBHs are modeled 
as sink particles that can accrete gas particles, allowing us 
to constrain the amount of mass accreted onto each MBH, 
and the orientation of the MBH spins relative to the angular 
momentum of the accreted gas. 

The paper is organized as follows: in Section 2 we focus 
on the SPH simulations, and describe the semi-analytical 
algorithm that evolves the MBH spins; in Section 3 we illus- 
trate our results on the MBH alignment and our prediction 
for the recoil velocity of the MBH remnant; in Section 4 we 
present our conclusions. 



2 NUMERICAL METHODS 
2.1 SPH simulations 

We follow the dynamics of MBH pairs in nuclear discs us- 
ing numerical simulations run with the N-Body/SPH code 
GADGET (Springel, Yoshida & White 2001), upgraded to 
include the accretion physics. The simulations discussed in 
this paper are the same as presented in Dotti et al. (2009). 
Here we give a short summary of the initial conditions for 
the different runs. For a more detailed discussion, we defer 
the reader to Dotti et al. (2009). 
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In our models, two MBHs are placed in the plane of a 
massive circumnuclear gaseous disc, embedded in a larger 
stellar spheroid. The gaseous disc is modeled with « 2 x 10 6 
particles, has a total mass Mm BC = 10 8 M Q , and follows a 
Mestel surface density profile S(i?) oc R~ , where R is the 
radial distance projected into the disc plane. The outer ra- 
dius of the disc is 100 pc. The massive disc is rotationally 
supported in R and has a vertical thickness of 8 pc. The 
internal energy per unit mass of the SPH particles scales as 
u(R) oc R~ 2 ^ A , where the value of the temperature at the 
outer radius of the disc has been set in order to have the 
Toomre parameter (Toomre 1964) Q > 3 everywhere, pre- 
venting the fragmentation of the disc (the average value of 
Q over the disc surface is « 10). Gas is evolved assuming a 
polytropic equation of state with index 7 = 5/3 or 7 = 7/5. 
In the former case, the runs are denoted by "H" and are 
termed "hot" as the temperature is proportional to a higher 
power of density than in the latter class of runs ("cold" 
cases, runs denoted by "C"). The cold case has been shown 
to provide a good approximation to a gas of solar metal- 
licity heated by a starburst (Spaans & Silk 2000; Klessen, 
Spaans, & Jappsen 2007). The hot case instead corresponds 
to an adiabatic monoatomic g if radiative cooling were 
completely suppressed during the merger, for example as a 
result of radiative heating after gas accretion onto the MBHs 
(Mayer et al. 2007). 

The spheroidal component (bulge) is modeled with 10 5 
collisionless particles, initially distributed as a Plummer 
sphere with a total mass MBui ge (= 6.98Moisc). The mass 
of the bulge within 100 pc is five times the mass of the disc, 
as suggested by Downes & Solomon (1998). 

The two MBHs (Mi and M 2 ) are equal in mass (M BH = 
4x 10 6 M ). The initial separation of the MBHs is 50 pc. Mi, 
called primary for reference, is placed at rest at the centre of 
the circumnuclear disc. M2, termed secondary, is moving on 
an initially eccentric (eo — 0.7) counterrotating (retrograde 
MBH, "R" runs) or corotating (prograde MBH, "P" runs) 
orbit with respect to the circumnuclear disc. Given the large 
masses of the disc and the bulge, the dynamics of the moving 
MBH (M2) is unaffected by the presence of Mi until the 
MBHs form a gravitationally bound system. 

We allow the gas particles to be accreted onto the MBHs 
if the following two criteria are fulfilled: 

• the sum of the kinetic and internal energy of the gas par- 
ticle is lower than 6-times the modulus of its gravitational 
energy (all the energies are computed with respect to each 
MBH); 

• the total mass accreted per unit time onto the MBH ev- 
ery timestep is lower than the accretion rate corresponding 
to the Eddington luminosity (Z/Edd) computed assuming a 
radiative efficiency (e) of 10%. 

The parameter 6 is a constant that defines the degree 
at which a particle is bound to the MBH in order to be ac- 
creted. We set b — 0.3. Note that due to the nature of the 
above criteria, the gas particles can accrete onto the MBHs 
only if the time- varying Bondi-Hoyle-Lyttleton radius is re- 
solved in the simulations. 

Each gas particle accreted by the MBH carries with 
it angular momentum. From the properties of the accreted 
particles we can compute, as a function of time, the mass 



Table 1. Run parameters 



run 


prograde ? 


eo 


7 


HP 


yes 


0.7 


5/3, "hot" 


HR 


no 


0.7 


5/3, "hot" 


CP 


yes 


0.7 


7/5, "cold" 


CR 


no 


0.7 


7/5, "cold" 



accretion rate and the versor l c dgc, that defines the direction 
of the total angular momentum of the accreted particles. 

This information can be gathered only by performing 
very high resolution simulations. The gravitational soften- 
ing of the MBHs is 0.1 pc. The gravitational softening of 
the gas particles is set to the same value, in order to prevent 
numerical errors. This is also the spatial resolution of the 
hydrodynamical force in the highest density regions^. The 
gravitational softening of the collisionless particles forming 
the bulge is 1 pc, in order to prevent two body interac- 
tions between gas particles and artificially massive stars. The 
main input parameters of our simulations are summarized 
in Table 1. 

2.2 Semi-analytical Bardeen-Petterson effect 

We use the MBH accretion histories obtained from our SPH 
simulations to follow the evolution of each MBH spin vector, 
Jbh = (aGM| H /c)JBH, where < a < 1 is the dimension- 
less spin parameter and Jbh is the spin versor. The scheme 
we adopt to study the spin evolution is based on the model 
recently developed by Perego et al. (2009). Here we summa- 
rize this algorithm. 

We assume that during any accretion event recorded in 
our SPH simulations, the inflowing gas forms a geometri- 
cally thin/optically thick a-disc (Shakura & Sunyaev 1973) 
on milli-parsec scales (not resolved in the simulation), and 
that the outer disc orientation is defined by the unit vec- 
tor ledge- The evolution of the a-disc is related to the radial 
viscosity v\ and the vertical viscosity v 2 '- v\ is the standard 
radial shear viscosity while ^2 is the vertical shear viscosity 
associated to the diffusion of vertical warps through the disc. 
The two viscosities can be described in terms of two different 
dimensionless viscosity parameters, a.\ and 02, through the 
relations ^1,2 = ai,2-Hc s , where H is the disc vertical scale 
height and c s is the sound speed of the gas in the accre- 
tion disc. We further assume «2 = /2/(2«i), with ot\ = 0.1 
and fi = 0.6 (Lodato & Pringle 2007). We assume power 
law profiles for the two viscosities, ^1,2 oc R :i ^ 4 , as in the 
Shakura & Sunyaev solution. 

As shown by Bardeen & Petterson (1975), if the or- 
bital angular momentum of the disc around the MBH is 
misaligned with respect to the MBH spin, the coupled action 
of viscosity and relativistic Lense-Thirring precession warps 

t The code computes the density of each SPH particle averaging 
over -/Vncigh = 32 neighbors. 
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the disc in its innermost region forcing the fluid to rotate 
in the equatorial plane of the spinning MBH. The timescale 
of propagation of the warp is short compared with the vis- 
cous/accretion timescale so that the deformed disc reaches 
an equilibrium profile that can be computed by solving the 
equation 



1 8 ( 3 dfi A , 
RdR{ UlER dR l ) + 



19/1 „ T d\ 
+ RdR\r* RL dR 



2G Jbh x L 
~ 2 



(1) 



where vr is the radial drift velocity, E is the surface density, 
and Q is the Keplerian velocity of the gas in the disc. L 
is the local angular momentum surface density of the disc, 
defined by its modulus L and the versor 1 that defines its 
direction. 

The boundary conditions to eq. 1 are the direction of L 
at the outer edge l c d go , the mass accretion rate (that fixes 
the magnitude of E), and the values of mass and spin of 
each MBH. All these values but the MBH spins are directly 
obtained from the SPH runs. In particular, the direction 
of the unit vector l e d g e is computed considering those SPH 
particles nearing the MBH gravitational sphere of influence 
that are accreted according to the criteria outlined in Section 
2.1. 

Also the MBH spin changes, not only because of accre- 
tion, but in response to its gravito-magnetic interaction with 
the disc on a timescale longer than the time scale of warp 
propagation (Perego et al. 2009). This interaction tends to 
reduce the degree of misalignment between the disc and the 
MBH spin, decreasing with time the angle between Jbh and 
lcdgo- The MBH spin evolution is followed by solving for the 
equation 

g(Jbh ,v,/n xi/n \ 47rG f L x Jbh 



dt 



MA(R lso )l(Riso) + 



4ttG f 

C 2 J A- 



R 2 



-dR. 



(2) 



The first term in eq. 2 accounts for the angular momentum 
deposited onto the MBH by the accreted particles at the in- 
nermost stable orbit (ISO), where A(i?rso) denotes the spe- 
cific angular momentum at 7?iso and l(iiiso) the unit vector 
parallel to Jbh, describing the warped disc according to the 
Bardeen-Petterson effect. The second term instead accounts 
for the gravo-magnetic interaction of the MBH spin with the 
warped disc. It modifies only the MBH spin direction (and 
not its modulus) , conserving the total angular momentum of 
the composite (MBH+disc) system (King et al. 2005). The 
integrand in eq. 2 peaks at the warp radius (i? W ar P ) where 
the disc deformation is the largest. t Eq. 2 incorporates two 
timescales: the accretion time related to the first right-hand 
term describing the e— folding increase of the spin modulus, 
and the shorter timescale of MBH spin alignment (Perego 
et al. 2009) 

-2/35 



Tal ~ 10° a 



5/7 



M BH 



4 x 10 6 M Q 



/ — 32/35 
Bdd y r i 



that will ensure a high degree of MBH-disc gravito-magnetic 
coupling during MBH inspiral, as we will show promptly in 
Section 3. In Eq 3 /Edd is the MBH luminosity in units of 

^Edd- 

We applied iteratively eq. 1 and 2 using inputs from the 
SPH simulation that give the values of the mass accretion 
rate, the MBH mass and the direction of l e d g e- The algorithm 
returns, as output, the spin vector, that is, its magnitude and 
direction. At each timestep our code therefore provides the 
angle between the spin vector of each MBH and the angular 
momentum vector of their relative orbit. 



3 RESULTS 

Figure 1 shows the time evolution of the relative angle 6 
between the spin of each MBH and the orbital angular mo- 



mentum of the MBH pair (L p 



iorblpair), for two selected 



(3) 



runs (CP and HR). The initial relative angle (ft) has been 
arbitrarily set to 2.5 radians (143°), while a has initially five 
different values (0.2, 0.4, 0.6, 0.8, and 1). 

There is a common trend in all the runs for both MBHs: 
MBHs with lower spins tend to align faster (as shown in 
Fig. 1 for t < 2 — 4 Myr) and are affected by changes in 
the plane of the accreting material to a larger extent (6 
changes rapidly with time and has more pronounced min- 
ima/maxima for lower a, see again Fig. 1). As indicated by 
eq. 3 a smaller spin modulus implies a shorter alignment 
time, and this explains the faster response of the MBH to 
orient its spin orthogonal to the plane of the accreted gas. 
A slowly-spinning MBH induces a weaker warp in the disc: 
the warp radius decreases with decreasing a and there the 
Lense-Thirring precession time is faster so that the MBH is 
more responsive to changes in the orientation of the accreted 
gas (see Perego et al. 2009 for details). 

The spin evolution depends also on the dynamical prop- 
erties of the MBHs and on the thermodynamics of the cir- 
cumnuclear disc. The effect of the initial orbital parameters 
is important during the first phase of orbital decay of the 
two MBHs, before they form a binary. We consider the two 
MBHs to be bound in a binary if the mass in gas and stars 
inside their orbits is lower than the mass of the binary. This 
happens when the separation between the two MBHs is ~ 5 
pc. The time at which the two MBHs form a binary (t hin ) 
is reported in Table 2 for each run. As described in detail 
in Section 3.1, M2 loses memory of its initial orbital pa- 
rameters before binding in a binary. As a consequence, the 
properties of accreting gas onto the MBHs after the for- 
mation of the binary are almost independent of the initial 
dynamical parameters of the pair. At this late stage of the 
orbital evolution, the gas accretion rate and the coherence of 
the accretion flows depend mostly on the thermodynamics 
of the circumnuclear disc. Summarizing, for t < ibin both 
dynamical and thermodynamical properties affect spin evo- 
lution, while for t > t hin the thermodynamical properties 
ultimately determine the final degree of spin alignment. 



t The exact definition of i? W arp is where the vertical viscous time 
R 2 1 in the disc is comparable to the Lense-Thirring precession 
time. Because i? W arp and the radius at which the disc is maximally 
deformed are comparable (Perego et al. 2009), wc simplify the 
notation in the paper using only -R war p- 



3.1 Effects of dynamics on spin alignment 

We note that for each MBH and every run, 8 initially 
(t < 4.5 Myr) decreases with time, but the alignment pro- 
cess is more efficient for Mi in both simulations. This delay 
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Figure 1. Upper panels: time evolution of the relative angle be- 
tween M2 spin and the orbital angular momentum of the MBH 
pair. Left (right) panel refers to runs CP (HR). The initial angle is 
arbitrarily set to 2.5 radians (close to anti-aligned) , and the initial 
spin parameter magnitudes varies between 0.2 (lighter colours) to 
1 (darker colours). Lower panels: same as upper panels for Mi. 



Table 2. Third column: MBH binary formation time. Fourth 
column: component parallel to L pa j r of the angular momentum 
of the gas particles accreted after the formation of the binary 
(AL Z ), normalized to its modulus (AL). Fifth column: average 
value of the angle between the MBH spins and L pa i r , after the 
formation of a binary. 



Run 


MBH 


tbin [Myr] 


AL Z /AL 


f (rad) 


CP 


Mi 


6.5 


>99.9% 


0.10 


CP 


Mi 


6.5 


>99.9% 


0.13 


CR 


Alt 


4.5 


>99.9% 


0.15 


CP, 


Ah 


4.5 


>99.9% 


0.16 


HP 


Mi 


7.5 


96.3% 


0.25 


HP 


Ah 


7.5 


94.9% 


0.23 


HR 


Alt 


4.5 


81.9% 


0.42 


HR 


Ah 


4.5 


77.9% 


0.32 



in the alignment of M2 is related to the orbital evolution of 
the orbiting MBH. In runs CP and HP, M2 is initially coro- 
tating with the circumnuclear disc on an eccentric orbit. 
Because of the eccentricity of the orbit, M2 has a non-zero 
relative velocity with respect to its local gas environment. As 
a consequence, the accretion rate onto M2 is initially lower 
than accretion rate onto Mi (Dotti et al. 2009). Dynamical 
friction exerted by the circumnuclear disc onto the orbiting 
MBH circularizes the orbit of M2 before the formation of 
a binary (Dotti, Colpi & Haardt 2006a; Dotti et al. 2007), 
so that the relative velocity between gas particles and M2 
decreases. After dynamical friction circularized the orbit of 
M2, the accretion rate onto M2 increases and becomes com- 



parable to the accretion rate onto Mi (Dotti et al. 2009). 
As a consequence the alignment of the spin of M2 becomes 
more efficient, and by the time a binary forms, 9 has similar 
values for Mi and M2 in the same run. 

For initially counterrotating MBHs (runs HR and CR), 
the effect of the dynamics onto the spin evolution of M2 
is more pronounced. Dynamical friction drags the orbiting 
MBH in the direction of the rotating gas, so that, before the 
formation of a binary, M2 starts to corotate with respect to 
the circumnuclear disc ("orbital angular momentum flip"; 
Dotti et al. 2009). In the counterrotating runs the ratio be- 
tween the accretion rate onto M2 before and after the angu- 
lar momentum flip can be < 0.15. As a consequence, during 
the first 2 — 3 Myrs, when the secondary moves on a retro- 
grade orbit, 9 does not change significantly (because of the 
low accretion rate), while it decreases efficiently only after 
the orbital angular momentum flip. 

3.2 Effects of gas thermodynamics on spin 
alignment 

Alignment occurs over a short time-scale, as indicated by 
the steep drop of 9 in Figure 1. Afterwards, 9 starts to oscil- 
late around an average value, different from run to run. Nu- 
merical noise due to the discrete nature of SPH calculations 
does not affect these oscillations. During each oscillation, 
the MBHs accrete tens of iVneigh. In particular, the average 
value of 9 and the amplitude of its oscillations are in general 
larger for hot runs (see Figure 1). We define 9{ as the angle 
between the MBH spins and L pa i r after the formation of a 
binary (9f = 9(t > tbin))- This new parameter is of key im- 
portance in the following discussion, since we assume that 
the distribution of 9f is representative of 9 at coalescence. 
The validity of this assumption is discussed in Section 3.3. 

The last column of Table 2 shows the average value 
of 9f of each MBH in every run. We note that 9{ is lower 
when the MBHs are embedded in colder discs. This is due to 
the properties of gas close to each MBH. For larger 7 (hot 
runs) the temperature of the gas in the overdense regions 
around each MBH is higher, and so is the pressure. As a 
consequence, the gas structures around each MBH (and the 
gas particles accreting onto the MBHs) are more pressure 
supported, spherical distributed, and with more isotropic 
velocities in runs HP and HR, while gas is more rotationally 
supported in runs CP and CR. This effect is quantified in 
Table 2. In the fourth column we report the component par- 
allel to Lpair of the angular momentum of the gas particles 
accreted after the formation of the binary (AL Z ), normal- 
ized to its modulus (AL). In cold runs, after the formation 
of the binary, streams of gas accreting onto the MBHs are 
extremely coherent (AL Z / AL > 99.9%). In hot runs the ac- 
creting particles have a larger degree of isotropy, resulting 
in less coherent accretion processes (AL Z /AL < 95%) and 
larger/more variable 9{. 

Since the time when a binary forms (ibin) is different for 
different runs, the time intervals (At) over which we average 
9{ are different. We decided to keep constant At for runs with 
the same polytropic index, but we used different At for cold 
and hot runs, in order to maximize the statistic. We chose 
At = 3.5 Myr for runs PC and RC, and At = 1 Myr for runs 
PH and RH. Averaging over different times does not affect 
the main results discussed above. As a check, we computed 
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Figure 2. Density distribution of pairs of the initial/final relative angles between MBH spins and orbital angular momentum. 

Left and right panels refer to MBHs embedded in cold and hot discs, respectively. Upper (lower) panels refer to the spin of M2 (Mi). 
Dark, medium, and light grey surfaces refer to high density regions encompassing 68.3%, 95.5%, and 99.7% of the realizations. 9\ has been 
sampled isotropically, and the dimcnsionlcss spin parameters (a) have been sampled from a constant probability distribution, between 
and 1. As discussed in the text, we average 8f over all the times after the formation of the MBH binary. 



9{ and At for the two MBHs in runs PC and RC averaging 
over only 1 Myr, and for every MBH/cold run combination 
we find t < 0.19 (10°) and AL Z /AL > 99%, consistent with 
the values reported in Table 2 for At = 3.5 Myr. 

We also note that physical processes not implemented in 
these simulations, such as star formation or feedback from 
supernovae, could decrease the degree of coherency of the 
accreting gas, possibly resulting in higher 8{. Furthermore, 
Lodato et al. 2009 have shown that star formation depletes 
the reservoir of gas in the vicinity of the MBHs, and can slow 



down the decay of the binary at sub-parsec separations. A 
detailed study of the interaction between star formation in 
the circumnuclear disc and the properties of the accreting 
gas is postponed to a future investigation. 

We estimated the efficiency of the alignment process 
over a large Monte Carlo sample of initial Jbh- For each 
MBH and each run, we selected 20,000 different initial val- 
ues of a, homogeneously distributed between and 1. For 
each value of a, we computed the three components of Jbh, 
assuming an initially isotropic distribution of the spins. We 
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evolved the initial condition for Jbh using the outputs of 
our simulations, as described in Section 2.2. As already dis- 
cussed in Section 3, the degree of alignment between MBH 
spins and L 

pair &t t ^> ibin 

is ultimately determined by the 
gas thermodynamics. As a consequence, we do not further 
analyse the dependence on the initial dynamics of M2, and 
focus mainly on the effect of the disc thermodynamics. The 
results from runs CP and CR have been combined in a single 
class (left panels in Figure 2). The same has been done for 
runs HP and HR (right panels in Figure 2). 

Figure 2 shows the density of realizations obtained with 
our statistical analysis, in the {Ovfit) plane. Dark, medium, 
and light grey surfaces refer to regions of decreasing density, 
encompassing 68.3%, 95.5%, and 99.7% of the realizations. 
We note that the alignment process is efficient independently 
of 9i. The lower density for 6i « and 6i « 7r is due to 
the initial isotropic distribution of the spins, and is totally 
unrelated to the alignment process. As already discussed 
above, alignment is more efficient for MBHs in cold discs. In 
these runs (left panels of Figure 2) 68.3% of the realizations 
have a final angle between the two MBHs and the orbital 
angular momentum of the pair 9t < 0.1(6°) while 68.3% 
of the realizations in runs HP and HR have 8{ < 0.5 (29°). 
There are a few % of the realizations with "large" final angles 
(0 f > 0.5 (29°)) in every run. 

3.3 Recoil distributions 

In this Section we assume that the two MBHs can reach 
coalescence, and we use the distributions of 9{ for Mi and 
M2 shown in Figure 2 in order to compute distributions 
of recoil velocities for the MBH remnant. We assume also 
that the distributions of 9{ we obtained are representative 
of the relative angle between the MBH spins and L pa i r dur- 
ing last phase of orbital decay, when the two MBHs lose 
efficiently orbital energy and angular momentum due to the 
gravitational wave emission. These assumptions are neces- 
sary because our simulations (spatial resolution w 0.1 pc) 
can not follow the evolution of the MBHs down to sepa- 
rations where gravitational waves dominate the dynamics. 
The two assumptions are valid if one of the following re- 
quirements is fulfilled: 

• The gas accretes onto the two MBHs in a coherent way un- 
til star formation and/or AGN feedback deplete the galactic 
nucleus of gas, and no further accretion events (i.e. due to 
tidal stripping of stars) change significantly the direction of 
the MBH spins; 

• The dynamical interaction between the binary and the 
gas creates a low density region (the so called "gap" , Gould 
& Rix 2000), reducing/halting accretion onto the MBHs 
(Milosavljevic & Phinney 2005; Dotti et al. 2006b; Hayasaki, 
Mineshige, Sudou 2007; Cuadra et al. 2009) so that the spins 
of the two MBHs do not change significantly when the bi- 
nary separation is < 0.1 pc; 

• After forming a binary, the two MBH can reach the final 
coalescence in a short time (< 10 Myr), so that further ac- 
cretion events do not have time to change significantly the 
MBH spin orientations. 

Numerical general relativistic computations show that 
the recoil velocity Vkick depends on the binary mass ratio 
q — M2/M1, on the dimensionless spin vectors of the pair 
ai and a2 (0 < a, < 1), and on the orbital parameters. This 
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Table 3. Recoil statistics. All velocities are in kms -1 



cold disc 


hot disc 


cold disc 


hot disc 


mean 


mean 


median 


median 



Fit CL 


40 


+50 

-40 


62 


+42 
-42 


23 


+ 13 

-13 


51 


+32 
-32 


Fit B 


39 


+45 
-39 


67 


+61 

-61 


24 


+ 12 

-12 


46 


+28 
-28 


Fit H 


33 


+27 
-27 


54 


+33 
-33 


23 


+7 
-7 


41 


+ 16 
-16 



information can be obtained from the analysis of our simula- 
tions. We use four different prescriptions from the literature 
to compute the recoil velocity of the MBH remnant, based 
on Campanelli et al. (2007) and Lousto & Zlochower (2009; 
fit CL), Baker et al. (2008; fit B), Herrmann et al. (2007; fit 
H), and Rezzolla et al. (2008; fit R). 

We use fit R as a consistency check, as this formula 
provides the recoil velocity for completely aligned configu- 
rations (6f from the simulations are close to, but not exactly 
zero), yielding lower limits for Vkick- When using fit R we 
adopt the spin magnitudes obtained from our simulations, 
further assuming that MBH spins are fully aligned with L pa i r 
at coalescence, and q = 1§. Expressions of the fitting formu- 
lae are detailed in the Appendix. 

The distributions of recoil velocities we obtain for cold 
(blue histograms) and hot (red histograms) discs are shown 
in Figure 3. For these three fitting formulae we report 
the distribution of recoil velocities we would obtain as- 
suming that the MBH spins are isotropically distributed 
(green lines). Because of the spin alignment discussed in Sec- 
tions 3.2 and 3.3, the recoil velocities we obtain analysing 
the results of our simulations are approximately one order 
of magnitude smaller than those predicted for isotropically 
distributed MBH spins, independently of the fitting formula 
we consider. Furthermore, the recoils obtained evolving the 
MBH pair in a cold circumnuclear disc are always a factor 
of « 2 smaller than the velocities obtained in the hot cases. 
This shift is due to the lower level of alignment between 
MBH spins and L pa i r for MBHs orbiting in hot discs, as dis- 
cussed in Section 3.2. The mean values of recoil velocities for 
these three fitting formulae and for different gas thermody- 
namics are shown the first two columns of Table 3. Because 
the mean values can be affected by the long tails of the recoil 
distributions at high velocities, we report also the median 
values in last two columns of the same table. 

Fit CL, fit B, and fit H give similar mean and median 
values, consistent within a factor of « 1- 1.3. The fraction 
of remnants with recoils larger than in cold (hot) runs with 
recoils larger than 400 kms -1 ("fast recoils") is 0.2% (8%) 
using fit CL. Fit B has the same fraction of fast recoils in 
cold runs, and a lower fraction (0.2%) in cold runs. Fit H 
does not have any realization with such high recoils. 

The black histogram in the upper panel shows the dis- 

§ The MBH mass ratio in our simulations is always between 0.9 
and 1. Assuming q = 1 does not affect our results. 
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Figure 3. Distribution of recoil velocities. The blue (red) histogram is computed from the distribution of spins we obtain from our 
simulations, after the alignment of the spins in a cold (hot) circumnuclear disc. The green curves refer to recoil velocities obtained 
assuming the spins of the two MBHs to be isotropically distributed. Upper, middle, and lower panels refer to the results obtained using 
fit H, fit CL, and fit B, respectively. In the upper panel the black histogram shows as a comparison the distribution of recoil velocity 
obtained using fit R, assuming complete alignment between the MBH spins and L pa i r . In this case both the results of cold and hot 
runs have been considered in a single histogram. In all the histograms the mass ratio between the MBHs (q) is obtained from our SPH 
simulations. 



tributions of recoil velocities obtained using fit R. In this 
case we considered both the results of cold and hot runs in 
a single histogram. Mean and median values for the recoils 
are w lOkms . Such low values follow from using the dis- 
tributions of a,i that we obtain from our simulations. After 
the formation of a binary, the MBHs in our runs have spin 
magnitudes 0.3 < a% < 0.9. If instead we assumed a homoge- 
neous distribution of spins between and 1, fit R would pre- 
dict a kick distribution with a peak at « 100 kms" 1 , a sharp 
cutoff at higher velocities, and a long tail at lower values. As 
expected, fit R gives a lower limit for the recoil velocities. 



The recoil distribution obtained with this last prescription 
peaks at velocities which are only m 2 times smaller than 
those where the cold runs peak, when using the other three 
fits. This confirms that these fitting formulae well describe 
quasi-aligned configurations. 



4 DISCUSSION 

In this paper, we traced for the first time the evolution of 
the spin vectors of MBHs orbiting inside a massive circum- 
nuclear gas disc. Our SPH simulations have sufficiently high 
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resolution to probe the hydrodynamics of fluid particles and 
the accretion physics near the gravitational sphere of influ- 
ence of the MBHs. An ad-hoc algorithm designed for track- 
ing the gravo-magnetic coupling between the MBH spin and 
the small-scale accretion disc is then implemented in the 
code. We find that: 

• When evolving a in dense, rotationally supported, struc- 
ture such as a circumnuclear disc, MBHs in a pair align their 
spins (Jbhi j2 ) to the pair orbital angular momentum (L pa ir) 
well before the two MBHs bind in a Keplerian binary, and 
independently of the MBHs initial orbital parameters. For a 
run with M2 initially on a retrograde orbit, the spin of the 
secondary aligns efficiently only after the "orbital angular 
momentum flip". 

• The average angle between Jbh! 2 an d L pa i r after the bi- 
nary formation (9f) depends on the thermodynamics of the 
massive circumnuclear discs. 6t is lower if the MBHs are em- 
bedded in colder discs (with a polytropic index 7 = 7/5), 
with respect to hotter discs (7 = 5/3); 

• After the formation of a binary, the two MBHs accrete gas 
with the same dynamical and thermodynamical properties. 
As a consequence, even the angle between the two small pro- 
jections of Jbh! 2 m the orbital plane decreases. This further 
reduces the recoil velocity of the MBH remnant. The degree 
of alignment between the two spins and between each spin 
and Lp a i r is preserved (or even increased) by spin-spin and 
spin-orbit interactions until the plunge phase (Schnittman 
2004; Herrmann et al. 2009); 

• Due to the efficient alignment between Jbh 1j2 an d Lpair, 
the expected recoil velocities (Vkick) at the MBH coalescence 
is, on average, one order of magnitude lower than those ex- 
pected for randomly oriented MBH spins. The thermody- 
namical properties of the environment affect the degree of 
alignment and, as a consequence, the expected recoil veloc- 
ities. Vkick is lower (by a factor of 1.5-2.2) for lower values 
of 7; 

• Assuming the same distribution of 9{, the recoil velocity 
distributions obtained using different prescriptions are very 
similar. The three fitting formulae used predict the same 
mean and median recoil velocities (< 70kms~ 1 ) within a 
factor of ~ 1 — 1.3, with less than few percent of realizations 
having Vkick > 400 km s -1 . 

The distributions of recoil velocities that we find have 
important consequences for retention of MBHs in galactic 
nuclei. When MBH binaries form and evolve in gas-rich ma- 
jor mergers, we predict the recoil velocity to be, on average, 
well below the escape speed from low-redshift galaxies. In- 
deed, because of the extreme efficiency of the spin align- 
ment process, the recoil velocities are likely unimportant 
even for high-z proto-galactic building blocks. Volonteri & 
Rees (2006) and Volonteri (2007) discussed how strong re- 
coils can affect the early growth of MBHs at the highest 
redshifts. In a forthcoming paper we will update our cal- 
culations and determine the impact of low recoils on the 
growth of MBHs in galaxies. 

Our simple treatment of thermodynamics and the ab- 
sence of any prescription for star-formation and supernovae 
feedback in our simulations could, in principle, overestimate 
the degree of coherency in the gas flows accreting onto 
the MBHs. Furthermore, our finite resolution prevent us 
to study the fragmentation of the accretion discs forming 
around the MBHs, that could result in a sequence of short 



and randomly oriented accretion events (King & Pringle 
2006). We will investigate interaction between star forma- 
tion in the circumnuclear disc and the properties of the ac- 
creting gas in a forthcoming study. 
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APPENDIX A: FITTING FORMULAE FOR 
RECOIL VELOCITIES 

Campanelli et al. (2007) and Lousto & Zlochower (2009; 
fit CL) propose the following fitting formula for the post- 
coalescence recoil of a MBH remnant: 



Vkick = ^v'L+v\+ 2v m v± cos(£) + ~v 2 
v m = Ar, 2 y/l-Ar 1 {l + Bri), 



v± = 



Hrf ( || || \ 

(iTg) r - qa v 

K v 2 
(1 + 9) 



cos(0 — Oo)|a^ — qa.2 



(4) 
(5) 
(6) 

(7) 



where A = 1.2 x 10 4 kms _1 , B = -0.93, H = 6900 kms" 1 , 
K — 6.0 x 10 4 kms _1 , 77 = g/(l + q) 2 is the symmetric 
mass ratio and £ measures the angle between the unequal 
mass and the spin contribution to the recoil velocity in the 
orbital plane. We assumed £ = 145°, as suggested by Lousto 
& Zlochower. The components of the spins of the two MBHs 
are: 



± 



0-2 



ai sin(6 l i) 

ai cos(#i) 

C12 sin(6 l 2) 

£J2 cos(#2), 



where the indices || and _L refer to projections parallel and 
perpendicular to the orbital angular momentum, respec- 
tively, and 81 (62) refers to 6{ for the primary (secondary) 
MBH. In eq. 7, O is the angle between (a^ — q&i) and the 
separation vector at coalescence, and 0o depends on the 
initial separation between the holes. Since Oo is unknown, 
for this exploration we assume a flat distribution of O — Oo 
between and 2tv. 
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Baker et al. (2008; fit B) propose instead the following 
fitting formula for v\\ : 

v\\ = (Y+qj cos ^ 1 ~ ~ 1 a 2 cos(0 2 - $2) ) ,(8) 

where cj>i (fo) is the angle between (a 2 ) and a fixed 
reference direction. Following Baker et al. (2008), $1 = 5>(g) 
and <E>2 = 3>(l/g). Because in our simulations q « 1, we 
fixed $1 = $2- We further assume a flat distribution of 
$1 between and 2tv. In this case, A = 1.35 x 10 4 kms _1 , 
B = -1.48, H = 7540 kms" 1 , and K = 2.4 x 10 5 kms _1 . 

Herrmann et al. (2007; fit H) formulate the recoil veloc- 
ity of a MBH remnant as a function of a different angle 9n , 
i.e. the angle between L pa i r and 

s -"'(t-Hr)' m 

where M = Mi + M2. Assuming that L pa i r is aligned with 
the z direction, they find that the Cartesian component of 
the recoil velocity follow: 

V x = C H x cos(6 u ), (10) 
V y = CoHyCosiOn), (11) 
V z = C K z sm{6u), (12) 

where Co = Eg 2 /(M 2 (1 + g) 4 ), and the best fitting parame- 
ters are H x = 2.1 xlO 3 H y = 7.3xlO a , and K z = 2.1xl0 4 T 
The last fitting formula we use to compute Vkick has 
been proposed by Rezzolla et al. (2008; fit R). In their study 
they consider only equal-mass MBHs with spins aligned 
with Lp a i r . They find: 

Vkick = \ci(ai - 02) + c 2 (a 2 - a|)|, (13) 

where ci = —220.97 and C2 = 45.52. Eq. 13 provides recoil 
velocities for completely aligned configurations. 
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